This guide is meant to assist the interpretation of analyses carried out within “Diebold and Gauthier et al. Effect of felzartamab anti-CD38 treatment on the molecular phenotype of antibody-mediated rejection in kidney transplant biopsies. under review. Nature Medicine”. This manusript is currently under peer-review and has not yet been accepted for publication. This guide is a companion to the GitHub repository (https://github.com/TSI-PTG/CD38-effect-of-treatment)
This analysis
assesses how felzartamab therapy affects molecular scores measured with
gene expression data from biopsies from kidney allografts taken from
patients with antibody-mediated rejection. The assumptions of normality
and heterogeneity of variance differ across molecular scores, thus we
avail to aligned-rank transformation to carry out a non-parametric
ANOVA. This is powerful in this context because it allows us to test for
the interactive effect of treatment (i.e., placebo vs felzartamab) and
time (i.e., change in scores from biopsy to biopsy). Aligned-rank
transformation also allows for us to specify a mixed-model to handle
repeated measurements (i.e., biopsies) across patients.
References:
• Wobbrock J, F.L., Gergle D, Higgins J The Aligned Rank Transform for Nonparametric Factorial Analyses Using Only ANOVA Procedures. In Proceedings of the ACM Conference on Human Factors in Computing Systems (CHI ’11) 143–146. doi:10.1145/1978942.1978963, https://depts.washington.edu/acelab/proj/art/.(2011).
• Kay M, E.L., Higgins J, Wobbrock J. ARTool: Aligned Rank Transform for Nonparametric Factorial ANOVAs. doi:10.5281/zenodo.594511, R package version 0.11.1, https://github.com/mjskay/ARTool. (2021).
• Elkin, L.A., Kay, M., Higgins, J.J. & Wobbrock, J.O. An Aligned Rank Transform Procedure for Multifactor Contrast Tests. in The 34th Annual ACM Symposium on User Interface Software and Technology 754–768 (Association for Computing Machinery, Virtual Event, USA, 2021).
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# Load necessary libraries
library(tidyverse)
## -- Attaching core tidyverse packages -------------------- tidyverse 2.0.0 --
## v dplyr 1.1.4 v readr 2.1.5
## v forcats 1.0.0 v stringr 1.5.1
## v ggplot2 3.5.1 v tibble 3.2.1
## v lubridate 1.9.3 v tidyr 1.3.1
## v purrr 1.0.2
## -- Conflicts -------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom)
library(rstatix)
##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filter
library(flextable)
##
## Attaching package: 'flextable'
##
## The following object is masked from 'package:purrr':
##
## compose
library(officer)
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(multcomp)
## Loading required package: mvtnorm
## Loading required package: survival
## Loading required package: TH.data
## Loading required package: MASS
##
## Attaching package: 'MASS'
##
## The following object is masked from 'package:rstatix':
##
## select
##
## The following object is masked from 'package:dplyr':
##
## select
##
##
## Attaching package: 'TH.data'
##
## The following object is masked from 'package:MASS':
##
## geyser
library(ARTool)
library(rcompanion)
library(Biobase)
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
##
## The following object is masked from 'package:flextable':
##
## width
##
## The following objects are masked from 'package:lubridate':
##
## intersect, setdiff, union
##
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
##
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
##
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter,
## Find, get, grep, grepl, intersect, is.unsorted, lapply, Map,
## mapply, match, mget, order, paste, pmax, pmax.int, pmin,
## pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
## setdiff, table, tapply, union, unique, unsplit, which.max,
## which.min
##
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
library(ggpubr)
##
## Attaching package: 'ggpubr'
##
## The following objects are masked from 'package:flextable':
##
## border, font, rotate
library(patchwork)
##
## Attaching package: 'patchwork'
##
## The following object is masked from 'package:MASS':
##
## area
library(ggprism)
# Custom operators
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Suppress dplyr reframe info
options(dplyr.reframe.inform = FALSE)
# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/results/flextable_artANOVA.RData")
# load plot data
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")
# Seed for reproducibility
set.seed(42)
Now we
define the categories of molecular scores for easier organization
downstream.
# Feature categories
vars_cfDNA <- c("cfDNA_cpml")
vars_abmr <- c("ABMRpm", "ggt0", "ptcgt0", "DSAST", "NKB")
vars_tcmr <- c("TCMRt", "tgt1", "igt1", "QCAT", "TCB")
vars_injury <- c("IRRAT30", "IRITD3", "IRITD5", "cigt1", "ctgt1")
vars_parenchyma <- c("KT1", "KT2")
vars_macrophage <- c("AMAT1", "QCMAT")
# DEFINE VARIABLES TO ASSESS ####
vars <- c(vars_cfDNA, vars_abmr, vars_tcmr, vars_macrophage, vars_injury, vars_parenchyma)
Now we
process the data so we can iteratively carry out the analyses on each
score.
# DEFINE THE DATA ####
data <- set %>%
pData() %>%
dplyr::select(Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, all_of(vars)) %>%
dplyr::filter(Patient %nin% c(15, 18)) %>%
left_join(., summarise(., sample_pairs = n(), .by = Felzartamab_Followup), by = "Felzartamab_Followup") %>%
relocate(sample_pairs, .after = Felzartamab_Followup) %>%
arrange(Felzartamab, Patient, Group)
# WRANGLE THE PHENOTYPE DATA ####
data_00 <- data %>%
expand_grid(category = c("cfDNA", "ABMR", "TCMR", "macrophage", "injury", "parenchyma")) %>%
nest(.by = category) %>%
mutate(
features = map(
category,
function(category) {
if (category == "cfDNA") {
vars_cfDNA
} else if (category == "ABMR") {
vars_abmr
} else if (category == "TCMR") {
vars_tcmr
} else if (category == "macrophage") {
vars_macrophage
} else if (category == "injury") {
vars_injury
} else if (category == "parenchyma") {
vars_parenchyma
}
}
),
data = pmap(
list(features, data),
function(features, data) {
data %>%
dplyr::select(
Center, Patient, Felzartamab, Group, Followup, Felzartamab_Group, Felzartamab_Followup, sample_pairs,
all_of(features)
)
}
)
)
# WRANGLE THE DATA FOR UNIVARIATE TESTS ####
data_01 <- data_00 %>%
mutate(
data_univariate = pmap(
list(data, features),
function(data, features) {
data %>%
pivot_longer(
cols = all_of(features),
names_to = "variable",
values_to = "value"
) %>%
nest(.by = variable) %>%
mutate(
variable = variable %>%
factor(
levels = vars
),
annotation = case_when(
variable %in% vars_cfDNA ~ "cfDNA",
variable %in% vars_tcmr ~ "TCMR-related",
variable %in% vars_abmr ~ "ABMR-related",
variable %in% vars_macrophage ~ "macrophage-related",
variable %in% vars_injury ~ "injury-related",
variable %in% vars_parenchyma ~ "parenchyma-related",
TRUE ~ " "
) %>%
factor(
levels = c(
"cfDNA",
"ABMR-related",
"TCMR-related",
"macrophage-related",
"injury-related",
"parenchyma-related",
"archetypes",
"rejection PC",
"injury PC"
)
),
score = case_when(
variable == "cfDNA_cpml" ~ "Donor-derived cell-free DNA (dd-cfDNA, cp/mL)",
variable == "TCMRt" ~ "TCMR classifier (TCMRProb)",
variable == "TCB" ~ "T cell burden (TCB)",
variable == "tgt1" ~ "Tubulitis classifier (t>1Prob)",
variable == "igt1" ~ "Interstitial infiltrate classifier (i>1Prob)",
variable == "QCAT" ~ "Cytotoxic T cell-associated (QCAT)",
variable == "ABMRpm" ~ "ABMR classifier (ABMRProb)",
variable == "DSAST" ~ "DSA-selective transcripts (DSAST)",
variable == "NKB" ~ "NK cell burden (NKB)",
variable == "ggt0" ~ "Glomerulitis classifier (g>0Prob)",
variable == "ptcgt0" ~ "Peritubular capillaritis classifier (ptc>0Prob)",
variable == "AMAT1" ~ "Alternatively activated macrophage (AMAT1)",
variable == "QCMAT" ~ "Constitutive macrophage (QCMAT)",
variable == "IRITD3" ~ "Injury-repair induced, day 3 (IRITD3)",
variable == "IRITD5" ~ "Injury-repair induced, day 5 (IRITD5)",
variable == "IRRAT30" ~ "Injury-repair associated (IRRAT30)",
variable == "cigt1" ~ "Fibrosis classifier (ci>1Prob)",
variable == "ctgt1" ~ "Atrophy classifier (ct>1Prob)",
variable == "KT1" ~ "Kidney parenchymal (KT1)",
variable == "KT2" ~ "Kidney parenchymal - no solute carriers (KT2)"
), .before = 1
) %>%
arrange(annotation, variable)
}
)
) %>%
dplyr::select(category, data_univariate) %>%
unnest(data_univariate) %>%
mutate(n_cat = score %>% unique() %>% length(), .by = category, .after = variable)
This is
what the data look like after we have organized them. Each row of the
dataframe contains the data for each score, the category it belongs to,
its category annotation, its long-format name, and the name of the
variable in the data.
data_01 %>% print(n="all")
## # A tibble: 20 x 6
## category annotation score variable n_cat data
## <chr> <fct> <chr> <fct> <int> <list>
## 1 cfDNA cfDNA Donor-derived cel~ cfDNA_c~ 1 <tibble>
## 2 ABMR ABMR-related ABMR classifier (~ ABMRpm 5 <tibble>
## 3 ABMR ABMR-related Glomerulitis clas~ ggt0 5 <tibble>
## 4 ABMR ABMR-related Peritubular capil~ ptcgt0 5 <tibble>
## 5 ABMR ABMR-related DSA-selective tra~ DSAST 5 <tibble>
## 6 ABMR ABMR-related NK cell burden (N~ NKB 5 <tibble>
## 7 TCMR TCMR-related TCMR classifier (~ TCMRt 5 <tibble>
## 8 TCMR TCMR-related Tubulitis classif~ tgt1 5 <tibble>
## 9 TCMR TCMR-related Interstitial infi~ igt1 5 <tibble>
## 10 TCMR TCMR-related Cytotoxic T cell-~ QCAT 5 <tibble>
## 11 TCMR TCMR-related T cell burden (TC~ TCB 5 <tibble>
## 12 macrophage macrophage-related Alternatively act~ AMAT1 2 <tibble>
## 13 macrophage macrophage-related Constitutive macr~ QCMAT 2 <tibble>
## 14 injury injury-related Injury-repair ass~ IRRAT30 5 <tibble>
## 15 injury injury-related Injury-repair ind~ IRITD3 5 <tibble>
## 16 injury injury-related Injury-repair ind~ IRITD5 5 <tibble>
## 17 injury injury-related Fibrosis classifi~ cigt1 5 <tibble>
## 18 injury injury-related Atrophy classifie~ ctgt1 5 <tibble>
## 19 parenchyma parenchyma-related Kidney parenchyma~ KT1 2 <tibble>
## 20 parenchyma parenchyma-related Kidney parenchyma~ KT2 2 <tibble>
At this
stage we can summarize the median values for each score by treatment and
follow-up.
(Note that the median scores are used for general interpretation of
treatment effects because aligned-ranks, the data being assessed for
inference regarding the effects of treatment, are not easily
interpretable.)
# UNIVARIATE MEDIANS ####
data_02 <- data_01 %>%
mutate(
medians = map(
data,
function(data) {
data %>%
reframe(
median = value %>% median(),
IQR = value %>% IQR(),
sample_pairs = n(),
.by = c(Followup, Felzartamab, Felzartamab_Followup)
) %>%
arrange(Felzartamab_Followup)
}
),
medians_delta = map(
medians,
function(medians) {
medians %>%
reframe(
median_delta = combn(median, 2, diff) %>% as.numeric(),
IQR_delta = combn(IQR, 2, mean) %>% as.numeric(),
sample_pairs = sample_pairs,
.by = c(Felzartamab)
) %>%
mutate(
Followup_pairwise = rep(c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52"), 2) %>%
factor(levels = c("Baseline - Week24", "Baseline - Week52", "Week24 - Week52")),
.before = sample_pairs
) %>%
mutate(
median_delta_delta = combn(median_delta, 2, diff) %>% as.numeric(),
IQR_delta_delta = combn(IQR_delta, 2, mean) %>% as.numeric(),
.by = c(Followup_pairwise),
.before = sample_pairs
)
}
)
)
Let’s
take a look at what the median summaries look like for the ABMRprob
score.
data_02 %>%
dplyr::filter(variable == "ABMRpm") %>%
pull(medians) %>%
pluck(1) %>%
flextable::flextable()
Followup | Felzartamab | Felzartamab_Followup | median | IQR | sample_pairs |
|---|---|---|---|---|---|
Baseline | Placebo | Baseline_Placebo | 0.6695 | 0.35525 | 10 |
Week24 | Placebo | Week24_Placebo | 0.7740 | 0.39025 | 10 |
Week52 | Placebo | Week52_Placebo | 0.5185 | 0.49825 | 10 |
Baseline | Felzartamab | Baseline_Felzartamab | 0.7360 | 0.30775 | 10 |
Week24 | Felzartamab | Week24_Felzartamab | 0.1675 | 0.35700 | 10 |
Week52 | Felzartamab | Week52_Felzartamab | 0.7025 | 0.53875 | 10 |
Now we
can look at the change in median ABMRprob scores in the placebo and
felzartamab arms, as well as the difference in those medians (i.e., ΔΔ
median - this is our summarized treatment effect).
data_02 %>%
dplyr::filter(variable == "ABMRpm") %>%
pull(medians_delta) %>%
pluck(1) %>%
flextable::flextable()
Felzartamab | median_delta | IQR_delta | Followup_pairwise | median_delta_delta | IQR_delta_delta | sample_pairs |
|---|---|---|---|---|---|---|
Placebo | 0.1045 | 0.372750 | Baseline - Week24 | -0.6730 | 0.3525625 | 10 |
Placebo | -0.1510 | 0.426750 | Baseline - Week52 | 0.1175 | 0.4250000 | 10 |
Placebo | -0.2555 | 0.444250 | Week24 - Week52 | 0.7905 | 0.4460625 | 10 |
Felzartamab | -0.5685 | 0.332375 | Baseline - Week24 | -0.6730 | 0.3525625 | 10 |
Felzartamab | -0.0335 | 0.423250 | Baseline - Week52 | 0.1175 | 0.4250000 | 10 |
Felzartamab | 0.5350 | 0.447875 | Week24 - Week52 | 0.7905 | 0.4460625 | 10 |
Now we
carry out the analysis. This includes testing the specific contrasts for
the interaction of treatment and time (i.e., our effect of
treatment).
# UNIVARIATE NONPARAMETRIC TESTS ####
artANOVA <- data_02 %>%
mutate(
art = map(data, art, formula = value ~ Followup * Felzartamab + (1 | Patient)),
art_aov = map(art, anova, type = "II"),
art_aov_tidy = map(art_aov, tidy) %>% suppressWarnings(),
art_aov_contrast = map(
art,
art.con,
formula = "Followup:Felzartamab", adjust = "fdr", method = "pairwise", interaction = TRUE, response = "art"
),
art_aov_contrast_tidy = map(art_aov_contrast, tidy),
art_aov_contrast_cld = map(
art_aov_contrast_tidy,
function(art_aov_contrast_tidy) {
art_aov_contrast_tidy %>%
as.data.frame() %>%
cldList(adj.p.value ~ Followup:Felzartamab, data = .)
}
)
)
# CREATE FLEXTABLE SUMMARY OF ART MODELS ####
title_art <- paste("Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients")
artANOVA %>%
dplyr::select(annotation, n_cat, score, art_aov) %>%
unnest(everything()) %>%
dplyr::rename(p.value = `Pr(>F)`) %>%
mutate(FDR = p.value %>% p.adjust(method = "fdr"), .by = annotation) %>%
dplyr::select(annotation, score, Term, `F`, p.value, FDR) %>%
flextable::flextable() %>%
flextable::add_header_row(values = rep(title_art, ncol_keys(.))) %>%
flextable::merge_h(part = "header") %>%
flextable::merge_v(j = 1:2) %>%
flextable::fontsize(size = 8, part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::bg(i = ~ FDR < 0.05 & Term == "Followup:Felzartamab", j = 3:6, bg = "#fbff00") %>%
flextable::colformat_double(j = 2:3, digits = 2) %>%
flextable::colformat_double(j = 4:ncol_keys(.), digits = 3) %>%
flextable::border_remove() %>%
flextable::bold(part = "header") %>%
flextable::padding(padding = 0, part = "all") %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::autofit()
Table i. Non-parametric ANOVA (ART) of molecular scores in biopsies from treated vs untreated patients | |||||
|---|---|---|---|---|---|
annotation | score | Term | F | p.value | FDR |
cfDNA | Donor-derived cell-free DNA (dd-cfDNA, cp/mL) | Followup | 4.418 | 0.019 | 0.029 |
Felzartamab | 1.598 | 0.222 | 0.222 | ||
Followup:Felzartamab | 10.396 | 0.000 | 0.001 | ||
ABMR-related | ABMR classifier (ABMRProb) | Followup | 6.659 | 0.003 | 0.007 |
Felzartamab | 1.620 | 0.219 | 0.235 | ||
Followup:Felzartamab | 9.605 | 0.000 | 0.001 | ||
Glomerulitis classifier (g>0Prob) | Followup | 11.768 | 0.000 | 0.001 | |
Felzartamab | 1.511 | 0.235 | 0.235 | ||
Followup:Felzartamab | 10.743 | 0.000 | 0.001 | ||
Peritubular capillaritis classifier (ptc>0Prob) | Followup | 14.188 | 0.000 | 0.000 | |
Felzartamab | 3.178 | 0.092 | 0.125 | ||
Followup:Felzartamab | 7.749 | 0.002 | 0.004 | ||
DSA-selective transcripts (DSAST) | Followup | 10.508 | 0.000 | 0.001 | |
Felzartamab | 1.747 | 0.203 | 0.234 | ||
Followup:Felzartamab | 4.179 | 0.023 | 0.035 | ||
NK cell burden (NKB) | Followup | 6.512 | 0.004 | 0.007 | |
Felzartamab | 1.956 | 0.179 | 0.224 | ||
Followup:Felzartamab | 5.114 | 0.011 | 0.018 | ||
TCMR-related | TCMR classifier (TCMRProb) | Followup | 3.555 | 0.039 | 0.165 |
Felzartamab | 1.417 | 0.249 | 0.468 | ||
Followup:Felzartamab | 1.124 | 0.336 | 0.504 | ||
Tubulitis classifier (t>1Prob) | Followup | 1.016 | 0.372 | 0.506 | |
Felzartamab | 1.853 | 0.190 | 0.408 | ||
Followup:Felzartamab | 0.065 | 0.937 | 0.968 | ||
Interstitial infiltrate classifier (i>1Prob) | Followup | 3.414 | 0.044 | 0.165 | |
Felzartamab | 0.997 | 0.331 | 0.504 | ||
Followup:Felzartamab | 0.032 | 0.968 | 0.968 | ||
Cytotoxic T cell-associated (QCAT) | Followup | 5.693 | 0.007 | 0.107 | |
Felzartamab | 4.118 | 0.057 | 0.172 | ||
Followup:Felzartamab | 0.927 | 0.405 | 0.506 | ||
T cell burden (TCB) | Followup | 4.360 | 0.020 | 0.151 | |
Felzartamab | 2.868 | 0.108 | 0.269 | ||
Followup:Felzartamab | 0.301 | 0.742 | 0.856 | ||
macrophage-related | Alternatively activated macrophage (AMAT1) | Followup | 1.549 | 0.226 | 0.304 |
Felzartamab | 1.047 | 0.320 | 0.320 | ||
Followup:Felzartamab | 2.406 | 0.105 | 0.209 | ||
Constitutive macrophage (QCMAT) | Followup | 3.378 | 0.045 | 0.209 | |
Felzartamab | 3.100 | 0.095 | 0.209 | ||
Followup:Felzartamab | 1.426 | 0.253 | 0.304 | ||
injury-related | Injury-repair associated (IRRAT30) | Followup | 0.055 | 0.946 | 0.946 |
Felzartamab | 3.842 | 0.066 | 0.223 | ||
Followup:Felzartamab | 1.915 | 0.162 | 0.347 | ||
Injury-repair induced, day 3 (IRITD3) | Followup | 0.174 | 0.841 | 0.946 | |
Felzartamab | 3.268 | 0.087 | 0.223 | ||
Followup:Felzartamab | 3.125 | 0.056 | 0.223 | ||
Injury-repair induced, day 5 (IRITD5) | Followup | 0.073 | 0.930 | 0.946 | |
Felzartamab | 1.088 | 0.311 | 0.583 | ||
Followup:Felzartamab | 2.586 | 0.089 | 0.223 | ||
Fibrosis classifier (ci>1Prob) | Followup | 0.360 | 0.700 | 0.876 | |
Felzartamab | 3.482 | 0.078 | 0.223 | ||
Followup:Felzartamab | 0.938 | 0.401 | 0.668 | ||
Atrophy classifier (ct>1Prob) | Followup | 0.568 | 0.572 | 0.780 | |
Felzartamab | 3.297 | 0.086 | 0.223 | ||
Followup:Felzartamab | 0.602 | 0.553 | 0.780 | ||
parenchyma-related | Kidney parenchymal (KT1) | Followup | 0.968 | 0.389 | 0.518 |
Felzartamab | 1.628 | 0.218 | 0.518 | ||
Followup:Felzartamab | 0.605 | 0.552 | 0.552 | ||
Kidney parenchymal - no solute carriers (KT2) | Followup | 0.859 | 0.432 | 0.518 | |
Felzartamab | 0.911 | 0.352 | 0.518 | ||
Followup:Felzartamab | 1.038 | 0.365 | 0.518 | ||
Here we
present the effect of treatment (ΔΔ) for each score, along with the
effect of time (Δ) for placebo and felzartmab arms.
flextable_pairwise
Table 1. Effect of felzartamab on molecular scores | |||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Annotation | Score | Baseline - Week 24 | Week 24 - Week 52 | Baseline - Week 52 | |||||||||
Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | Δ Placebo | Δ Felzartamab | ΔΔ | ΔΔ FDR | ||
cfDNA | Donor-derived cell-free DNA (dd-cfDNA, cp/mL) | 0.50 (76.12) | -55.00 (30) | -55.50 (53.06) | 2e-04 | -3.50 (69.12) | 14.50 (25.5) | 18.00 (47.31) | 0.027 | -3.00 (51) | -40.50 (45.75) | -37.50 (48.38) | 0.045 |
ABMR-related | ABMR classifier (ABMRProb) | 0.10 (0.37) | -0.57 (0.33) | -0.67 (0.35) | 1e-03 | -0.26 (0.44) | 0.54 (0.45) | 0.79 (0.45) | 1e-03 | -0.15 (0.43) | -0.03 (0.42) | 0.12 (0.43) | 0.709 |
Glomerulitis classifier (g>0Prob) | -0.02 (0.23) | -0.31 (0.35) | -0.29 (0.29) | 3e-03 | -0.19 (0.34) | 0.31 (0.41) | 0.50 (0.37) | 2e-04 | -0.22 (0.3) | 0.00 (0.38) | 0.21 (0.34) | 0.277 | |
Peritubular capillaritis classifier (ptc>0Prob) | -0.12 (0.25) | -0.55 (0.27) | -0.43 (0.26) | 3e-03 | -0.10 (0.37) | 0.37 (0.31) | 0.47 (0.34) | 3e-03 | -0.22 (0.3) | -0.18 (0.35) | 0.04 (0.33) | 0.977 | |
DSA-selective transcripts (DSAST) | -0.09 (0.22) | -0.37 (0.24) | -0.28 (0.23) | 0.026 | -0.16 (0.32) | 0.16 (0.33) | 0.32 (0.33) | 0.026 | -0.24 (0.33) | -0.20 (0.25) | 0.04 (0.29) | 0.988 | |
NK cell burden (NKB) | 0.02 (0.22) | -0.78 (0.49) | -0.80 (0.36) | 0.014 | -0.19 (0.4) | 0.68 (0.59) | 0.87 (0.49) | 0.014 | -0.17 (0.46) | -0.10 (0.62) | 0.07 (0.54) | 0.972 | |
TCMR-related | TCMR classifier (TCMRProb) | -0.01 (0.02) | 0.00 (0) | 0.01 (0.01) | 0.358 | 0.00 (0.02) | 0.00 (0) | 0.00 (0.01) | 0.856 | -0.01 (0.02) | 0.00 (0) | 0.01 (0.01) | 0.358 |
Tubulitis classifier (t>1Prob) | -0.01 (0.06) | -0.02 (0.04) | -0.01 (0.05) | 1.000 | 0.00 (0.04) | 0.01 (0.02) | 0.01 (0.03) | 1.000 | -0.01 (0.05) | -0.01 (0.03) | -0.01 (0.04) | 1.000 | |
Interstitial infiltrate classifier (i>1Prob) | -0.01 (0.04) | -0.02 (0.03) | -0.01 (0.04) | 0.964 | 0.01 (0.04) | 0.01 (0.04) | 0.00 (0.04) | 0.964 | 0.00 (0.05) | -0.01 (0.06) | -0.01 (0.05) | 0.964 | |
Cytotoxic T cell-associated (QCAT) | -0.20 (0.41) | -0.70 (0.46) | -0.50 (0.43) | 0.425 | 0.10 (0.33) | 0.54 (0.57) | 0.44 (0.45) | 0.425 | -0.10 (0.46) | -0.16 (0.41) | -0.06 (0.44) | 0.871 | |
T cell burden (TCB) | -0.44 (0.71) | -0.79 (0.56) | -0.36 (0.64) | 0.825 | 0.28 (0.57) | 0.34 (0.6) | 0.06 (0.58) | 0.905 | -0.16 (0.88) | -0.45 (0.43) | -0.29 (0.65) | 0.825 | |
macrophage-related | Alternatively activated macrophage (AMAT1) | 0.08 (0.2) | -0.27 (0.45) | -0.35 (0.33) | 0.128 | 0.09 (0.23) | 0.11 (0.47) | 0.02 (0.35) | 0.812 | 0.18 (0.19) | -0.16 (0.32) | -0.33 (0.25) | 0.128 |
Constitutive macrophage (QCMAT) | 0.00 (0.15) | -0.24 (0.27) | -0.24 (0.21) | 0.340 | 0.01 (0.12) | 0.09 (0.27) | 0.08 (0.2) | 0.685 | 0.01 (0.11) | -0.15 (0.24) | -0.16 (0.18) | 0.348 | |
injury-related | Injury-repair associated (IRRAT30) | 0.18 (0.52) | 0.05 (0.4) | -0.13 (0.46) | 0.353 | 0.05 (0.66) | -0.11 (0.72) | -0.16 (0.69) | 0.353 | 0.23 (0.6) | -0.06 (0.48) | -0.29 (0.54) | 0.175 |
Injury-repair induced, day 3 (IRITD3) | 0.03 (0.19) | 0.04 (0.12) | 0.01 (0.15) | 0.945 | 0.03 (0.18) | -0.18 (0.18) | -0.20 (0.18) | 0.060 | 0.06 (0.19) | -0.13 (0.14) | -0.19 (0.16) | 0.060 | |
Injury-repair induced, day 5 (IRITD5) | -0.01 (0.17) | 0.04 (0.14) | 0.05 (0.15) | 0.932 | 0.07 (0.2) | -0.17 (0.21) | -0.24 (0.21) | 0.093 | 0.06 (0.18) | -0.13 (0.17) | -0.18 (0.18) | 0.093 | |
Fibrosis classifier (ci>1Prob) | -0.02 (0.29) | 0.13 (0.32) | 0.15 (0.3) | 0.434 | 0.06 (0.25) | -0.08 (0.53) | -0.14 (0.39) | 0.844 | 0.04 (0.25) | 0.05 (0.37) | 0.01 (0.31) | 0.434 | |
Atrophy classifier (ct>1Prob) | 0.01 (0.33) | 0.12 (0.31) | 0.12 (0.32) | 0.587 | 0.09 (0.27) | -0.06 (0.44) | -0.15 (0.36) | 0.587 | 0.10 (0.32) | 0.07 (0.41) | -0.03 (0.36) | 0.587 | |
parenchyma-related | Kidney parenchymal (KT1) | -0.05 (0.22) | -0.15 (0.23) | -0.09 (0.22) | 0.731 | -0.01 (0.11) | 0.10 (0.22) | 0.11 (0.17) | 0.705 | -0.06 (0.25) | -0.04 (0.18) | 0.02 (0.22) | 0.705 |
Kidney parenchymal - no solute carriers (KT2) | -0.09 (0.39) | -0.16 (0.37) | -0.07 (0.38) | 0.577 | -0.09 (0.2) | 0.14 (0.38) | 0.23 (0.29) | 0.484 | -0.18 (0.42) | -0.02 (0.31) | 0.16 (0.37) | 0.577 | |
Grey shading denotes ANOVA interactive effect FDR < 0.05 | |||||||||||||
# EXTRACT LEGEND FOR PLOTS ####
panel_legend <- plots_artANOVA %>%
dplyr::filter(variable == "AMAT1") %>%
pull(plot_violin) %>%
ggpubr::get_legend() %>%
ggpubr::as_ggplot() +
theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))
# MOLECULAR ABMR PANELS ####
panel_violin_abmr <- plots_artANOVA %>%
dplyr::filter(category %in% c("ABMR")) %>%
pull(plot_violin) %>%
wrap_plots(nrow = 1, ncol = 5) &
theme(
axis.text = element_text(size = 10, colour = "black"),
legend.position = "none",
plot.background = element_rect(fill = "grey95", colour = " white")
)
panel_violin_abmr <- panel_violin_abmr %>%
ggarrange(
labels = "B",
font.label = list(size = 25, face = "bold"),
legend = "none"
) %>%
ggpubr::annotate_figure(
top = text_grob("Effect of felzartamab on molecular ABMR activity scores",
face = "bold.italic",
size = 25,
hjust = 1.27
)
)
# MOLECULAR TCMR PANELS ####
panel_violin_tcmr <- plots_artANOVA %>%
dplyr::filter(category %in% c("TCMR")) %>%
pull(plot_violin) %>%
wrap_plots(nrow = 1, ncol = 5) &
theme(
axis.text = element_text(size = 10, colour = "black"),
legend.position = "none",
)
panel_violin_tcmr <- panel_violin_tcmr %>%
ggarrange(
labels = "C",
font.label = list(size = 25, face = "bold"),
legend = "none"
) %>%
ggpubr::annotate_figure(
top = text_grob("Effect of felzartamab on molecular TCMR activity scores",
face = "bold.italic",
size = 25,
hjust = 1.27
)
)
# COMBINED VIOLIN PANELS ####
panels_violin <- ggarrange(
panel_violin_abmr,
panel_violin_tcmr,
nrow = 2,
heights = c(1, 1)
)
ggarrange(
panel_legend,
panels_violin,
nrow = 2,
heights = c(0.125, 1)
)
This analysis
assesses the slopes for injury scores over time between placebo and
felzartamab arms. We used a linear mixed-model to account for repeated
measurements over time.
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# CRAN packages
library(tidyverse)
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
##
## lmer
## The following object is masked from 'package:stats':
##
## step
library(performance)
library(ggeffects)
library(flextable)
library(ggpubr)
library(patchwork)
library(ggprism)
# Bioconductor libraries
library(Biobase) # BiocManager::install("Biobase")
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
source("C:/R/CD38-effect-of-treatment/natmed/code/functions/get_slope_function.r")
# load reference set
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
# load refression plots
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_regression.RData")
# load violin plots
load("C:/R/CD38-effect-of-treatment/natmed/results/plots_artANOVA.RData")
# Seed for reproducibility
set.seed(42)
# DEFINE THE DATA ####
data <- set %>%
pData() %>%
dplyr::select(Center, Patient, Felzartamab, Group, IRRAT30, IRITD3, IRITD5) %>%
dplyr::filter(Patient %nin% c(15, 18)) %>%
mutate(
time = case_when(
Group == "Index" ~ 0,
Group == "FU1" ~ 0.46,
Group == "FU2" ~ 0.99
),
linetype = "solid" %>% factor()
)
Now we
fit the linear mixed-models.
# FIT MODELS ####
fit_IRRAT30 <- lmer(IRRAT30 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD3 <- lmer(IRITD3 ~ Felzartamab * time + (time | Patient), data)
fit_IRITD5 <- lmer(IRITD5 ~ Felzartamab * time + (time | Patient), data)
We check
the validatity of the model assumptions for IRRAT30.
# test model assumptions
performance::check_model(fit_IRRAT30)
performance::check_normality(fit_IRRAT30)
## OK: residuals appear as normally distributed (p = 0.605).
We check
the validatity of the model assumptions for IRITD3.
# test model assumptions
performance::check_model(fit_IRITD3)
performance::check_normality(fit_IRITD3)
## OK: residuals appear as normally distributed (p = 0.070).
We check
the validatity of the model assumptions for IRITD5.
# test model assumptions
performance::check_model(fit_IRITD5)
performance::check_normality(fit_IRITD5)
## OK: residuals appear as normally distributed (p = 0.090).
The next
step is to extract and summarize the model results.
# EXTRACT THE SLOPES FROM THE MODEL FITS ####
# IRRAT30
# Create list object into which to place model results
slopes_IRRAT30 <- list()
# Acute slopes
slopes_IRRAT30[[1]] <- get_slope(
.model_obj = fit_IRRAT30, .name = "Placebo",
.contrasts = c(0, 0, 1, 0) #
)
slopes_IRRAT30[[2]] <- get_slope(
.model_obj = fit_IRRAT30, .name = "Felzartamab",
.contrasts = c(0, 0, 1, 1)
)
slopes_IRRAT30[[3]] <- get_slope(
.model_obj = fit_IRRAT30, .name = "Intergroup Difference",
.contrasts = c(0, 0, 0, 1)
)
model_IRRAT30 <- slopes_IRRAT30 %>%
bind_rows() %>%
mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
dplyr::select(name, result, p_value, p_adjusted)
# IRITD3
# Create list object into which to place model results
slopes_IRITD3 <- list()
# Acute slopes
slopes_IRITD3[[1]] <- get_slope(
.model_obj = fit_IRITD3, .name = "Placebo",
.contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD3[[2]] <- get_slope(
.model_obj = fit_IRITD3, .name = "Felzartamab",
.contrasts = c(0, 0, 1, 1)
)
slopes_IRITD3[[3]] <- get_slope(
.model_obj = fit_IRITD3, .name = "Intergroup Difference",
.contrasts = c(0, 0, 0, 1)
)
model_IRITD3 <- slopes_IRITD3 %>%
bind_rows() %>%
mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
dplyr::select(name, result, p_value, p_adjusted)
# IRITD5
# Create list object into which to place model results
slopes_IRITD5 <- list()
# Acute slopes
slopes_IRITD5[[1]] <- get_slope(
.model_obj = fit_IRITD5, .name = "Placebo",
.contrasts = c(0, 0, 1, 0) #
)
slopes_IRITD5[[2]] <- get_slope(
.model_obj = fit_IRITD5, .name = "Felzartamab",
.contrasts = c(0, 0, 1, 1)
)
slopes_IRITD5[[3]] <- get_slope(
.model_obj = fit_IRITD5, .name = "Intergroup Difference",
.contrasts = c(0, 0, 0, 1)
)
model_IRITD5 <- slopes_IRITD5 %>%
bind_rows() %>%
mutate(p_adjusted = p.adjust(p_value, method = "fdr")) %>%
dplyr::select(name, result, p_value, p_adjusted)
IRRAT30
size_font <- 10
model_IRRAT30 %>%
dplyr::select(-p_value) %>%
flextable::qflextable() %>%
flextable::set_table_properties(layout = "autofit") %>%
flextable::set_header_labels(
name = "Group", result = "IRRAT30 Slope (95%CI)",
"p_adjusted" = "FDR"
) %>%
flextable::bold(i = NULL, part = "header") %>%
flextable::fontsize(size = size_font, part = "all")
Group | IRRAT30 Slope (95%CI) | FDR |
|---|---|---|
Placebo | 0.32 (-0.04 to 0.67) | 0.1133 |
Felzartamab | -0.29 (-0.64 to 0.07) | 0.1133 |
Intergroup Difference | -0.60 (-1.10 to -0.10) | 0.0546 |
IRITD3
model_IRITD3 %>%
dplyr::select(-p_value) %>%
flextable::qflextable() %>%
flextable::set_table_properties(layout = "autofit") %>%
flextable::set_header_labels(
name = "Group", result = "IRITD3 Slope (95%CI)",
"p_adjusted" = "FDR"
) %>%
flextable::bold(i = NULL, part = "header") %>%
flextable::fontsize(size = size_font, part = "all")
Group | IRITD3 Slope (95%CI) | FDR |
|---|---|---|
Placebo | 0.06 (-0.06 to 0.18) | 0.3453 |
Felzartamab | -0.14 (-0.26 to -0.02) | 0.0408 |
Intergroup Difference | -0.20 (-0.37 to -0.02) | 0.0408 |
IRITD5
model_IRITD5 %>%
dplyr::select(-p_value) %>%
flextable::qflextable() %>%
flextable::set_table_properties(layout = "autofit") %>%
flextable::set_header_labels(
name = "Group", result = "IRITD5 Slope (95%CI)",
"p_adjusted" = "FDR"
) %>%
flextable::bold(i = NULL, part = "header") %>%
flextable::fontsize(size = size_font, part = "all")
Group | IRITD5 Slope (95%CI) | FDR |
|---|---|---|
Placebo | 0.09 (-0.04 to 0.22) | 0.15720 |
Felzartamab | -0.13 (-0.26 to -0.01) | 0.05895 |
Intergroup Difference | -0.23 (-0.41 to -0.05) | 0.04200 |
# DEFINE LEGEND FOR VIOLIN PLOTS ####
panel_legend_violin <- plots_artANOVA %>%
dplyr::filter(variable == "AMAT1") %>%
pull(plot_violin) %>%
ggpubr::get_legend() %>%
ggpubr::as_ggplot() +
theme(plot.margin = unit(c(0, 0, -1, 0), "cm"))
# DEFINE LEGEND FOR REGRESSION PLOTS ####
panel_legend_regression <- plots_regression %>%
dplyr::filter(variable == "IRRAT30") %>%
pull(plot) %>%
ggpubr::get_legend() %>%
ggpubr::as_ggplot()
# DEFIN JOINT LEGEND FOR REGRESSION AND VIOLIN PLOTS ####
panel_legend_all_injury <- panel_legend_violin + panel_legend_regression
# DEFINE INJURY VIOLIN PLOTS ####
plots_violin <- plots_artANOVA %>%
dplyr::filter(category %in% c("injury")) %>%
pull(plot_violin)
# MAKE PANEL OF SLOPE AND VIOLIN PLOTS ####
plot_panels_all_injury <- patchwork::wrap_plots(
plots_violin[[1]], plots_violin[[2]], plots_violin[[3]],
plots_regression$plot[[1]], plots_regression$plot[[2]], plots_regression$plot[[3]],
plots_regression$grob_table[[1]], plots_regression$grob_table[[2]], plots_regression$grob_table[[3]],
nrow = 3,
ncol = 3,
guides = "collect",
heights = c(1, 1, 0.5)
) +
patchwork::plot_annotation(tag_levels = list(c(LETTERS[1:6], rep("", 3)))) &
theme(
legend.position = "none",
# plot.title = element_blank(),
axis.text = element_text(size = 10, colour = "black"),
axis.title = element_text(size = 12, colour = "black"),
plot.tag = element_text(size = 20, face = "bold", vjust = 1)
)
ggarrange(
panel_legend_all_injury,
plot_panels_all_injury,
nrow = 2,
heights = c(0.15, 1.5)
)
This analysis
assesses how felzartamab therapy affects gene expression measured in
kidney allograft biopsies taken from patients with antibody mediated
rejection. We use the Linear Models for Microarray Data (limma) pipeline
for applying empirical Bayes mediated t-tests on individual genes.
Specific contrasts to assess the interaction between treatment and time
are specified to estimate the effect of treatment on each gene.
References:
• Ritchie, M.E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse)
library(flextable)
library(officer)
# Bioconductor libraries
library(Biobase)
library(limma)
##
## Attaching package: 'limma'
## The following object is masked from 'package:BiocGenerics':
##
## plotMA
library(genefilter)
##
## Attaching package: 'genefilter'
## The following object is masked from 'package:patchwork':
##
## area
## The following object is masked from 'package:MASS':
##
## area
## The following object is masked from 'package:rstatix':
##
## Anova
## The following object is masked from 'package:readr':
##
## spec
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# Load data
load("C:/R/CD38-effect-of-treatment/natmed/data/cd38_3Sept24.RData")
load("C:/R/CD38-effect-of-treatment/natmed/data/affymap219.RData")
# Seed for reproducibility
set.seed(42)
# DEFINE THE SET ####
set00 <- set[, set$Patient %nin% c(15, 18)]
• We
filter the gene expression data by to remove genes with minimal variance
across the population (IQR filtering).
# IQR FILTER THE DATA ####
f1 <- function(x) (IQR(x) > 0.5)
ff <- filterfun(f1)
if (!exists("selected")) {
selected <- genefilter(set00, ff)
}
set01 <- set00[selected, ]
• We then
retain a single probeset for each gene based on whichever has the
highest mean expression in the population.
# KEEP UNIQUE GENES (keep probe with highest mean expression) ####
mean_exprs_by_probe <- set01 %>%
exprs() %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
mutate(mean_exprs = set01 %>%
exprs() %>% rowMeans(), .after = Symb)
genes <- mean_exprs_by_probe %>%
group_by(Symb) %>%
dplyr::slice_max(mean_exprs) %>%
dplyr::filter(Symb != "") %>%
distinct(Symb, .keep_all = TRUE) %>%
pull(AffyID)
set <- set01[featureNames(set01) %in% genes, ]
The
felzartamab trial was a phase II randomized control trial (RCT). While
RCTs are the ‘gold standard’ for accounting for confounders while
establishing causality (i.e., the effect of treatment), they are not
perfect, particularly when the number of study participants is low.
Thus, we pre-screened gene expression data to remove any genes that
differed between placebo and felzartamab arms at baseline. This was done
to reduce the potential for interpreting the effect of treatment on
genes with dissimilar baseline profiles.
# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set$Felzartamab_Followup %>% droplevels()
# BLOCK DESIGN week24 - baseline ####
design_block01 <- model.matrix(~ 0 + Felzartamab_Followup)
contrast_block_01 <- makeContrasts(
"x = (Felzartamab_FollowupBaseline_Felzartamab) -(Felzartamab_FollowupBaseline_Placebo)",
levels = design_block01
)
# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_block_1 <- limma::lmFit(set, design_block01)
cfit_block_1 <- limma::contrasts.fit(fit_block_1, contrast_block_01)
ebayes_block_1 <- limma::eBayes(cfit_block_1)
tab_block_1 <- limma::topTable(ebayes_block_1, adjust = "fdr", sort.by = "p", number = "all")
# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means_baseline <- fit_block_1 %>%
avearrays() %>%
data.frame() %>%
rownames_to_column("AffyID") %>%
tibble() %>%
mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
dplyr::select(AffyID, contains("Baseline"))
The top
20 genes that differ at baseline (sorted by p-value) between placebo and
felzartamab arms are summarized below.
(Note that all FDR are > 0.05. This is a consequence of differential expression across a large selection of genes (5,267 in this case) with small samples sizes (i.e., n = 10 by treatment group). Nonetheless, genes with p < 0.05 will be excluded from subsequent analyses).
# FORMAT TOPTABLES ####
baseline_deg <- tab_block_1 %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID") %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means_baseline, by = "AffyID") %>%
dplyr::select(
AffyID, Symb,
Baseline_Placebo, Baseline_Felzartamab,
t, logFC, P.Value, adj.P.Val,
) %>%
dplyr::rename(p = P.Value, FDR = adj.P.Val)
baseline_deg %>%
dplyr::slice_min(p, n = 20) %>%
mutate_at(vars(t, logFC, FDR), ~ formatC(.x, format = "f", digits = 2)) %>%
mutate_at(vars(p), ~ formatC(.x, format = "f", digits = 5)) %>%
flextable::flextable()
AffyID | Symb | Baseline_Placebo | Baseline_Felzartamab | t | logFC | p | FDR |
|---|---|---|---|---|---|---|---|
11727331_at | SPAG4 | 90 | 50 | -4.23 | -0.85 | 0.00008 | 0.33 |
11763725_a_at | HNRNPU | 366 | 665 | 4.10 | 0.86 | 0.00013 | 0.33 |
11718286_a_at | UBL3 | 74 | 142 | 3.66 | 0.94 | 0.00053 | 0.45 |
11732590_x_at | TRIM5 | 42 | 63 | 3.65 | 0.59 | 0.00055 | 0.45 |
11725671_a_at | IPO7 | 42 | 66 | 3.54 | 0.64 | 0.00079 | 0.45 |
11741431_a_at | KCNJ16 | 203 | 311 | 3.49 | 0.61 | 0.00091 | 0.45 |
11726915_at | SLCO4C1 | 87 | 176 | 3.48 | 1.02 | 0.00096 | 0.45 |
11723109_a_at | GPBP1 | 153 | 263 | 3.46 | 0.78 | 0.00099 | 0.45 |
11739744_x_at | MXI1 | 123 | 191 | 3.45 | 0.64 | 0.00102 | 0.45 |
11720890_a_at | AGL | 39 | 63 | 3.40 | 0.70 | 0.00121 | 0.45 |
11746431_a_at | FBXO11 | 28 | 45 | 3.38 | 0.68 | 0.00127 | 0.45 |
11751841_a_at | TNFRSF17 | 56 | 20 | -3.36 | -1.51 | 0.00134 | 0.45 |
11744367_a_at | ZMYND11 | 30 | 49 | 3.33 | 0.68 | 0.00150 | 0.45 |
11748857_a_at | AK3 | 371 | 578 | 3.32 | 0.64 | 0.00155 | 0.45 |
11727307_x_at | ZNF320 | 46 | 66 | 3.31 | 0.51 | 0.00160 | 0.45 |
11747648_a_at | TCAIM | 96 | 152 | 3.29 | 0.66 | 0.00168 | 0.45 |
11728966_s_at | ZNF28 | 29 | 67 | 3.27 | 1.21 | 0.00181 | 0.45 |
11723387_a_at | PTAR1 | 149 | 261 | 3.25 | 0.81 | 0.00187 | 0.45 |
11739771_a_at | OXR1 | 31 | 51 | 3.24 | 0.73 | 0.00194 | 0.45 |
11730038_at | GNG13 | 44 | 28 | -3.24 | -0.63 | 0.00196 | 0.45 |
We now
carry out differential expression analysis to estimate the treatment
effect on individual genes.
# DEFINE GENES SIMILAR AT BASELINE ####
genes_baseline <- baseline_deg %>%
dplyr::filter(p > 0.05) %>%
pull(AffyID)
# WRANGLE THE MEAN EXPRESSION DATA ####
mean_exprs_by_probe <- set01 %>%
exprs() %>%
as_tibble(rownames = "AffyID") %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb) %>% tibble(), ., by = "AffyID") %>%
mutate(
mean_exprs = set01 %>%
exprs() %>%
rowMeans(),
.after = Symb
)
genes <- mean_exprs_by_probe %>%
group_by(Symb) %>%
dplyr::slice_max(mean_exprs) %>%
dplyr::filter(Symb != "", AffyID %in% genes_baseline) %>%
distinct(Symb, .keep_all = TRUE) %>%
pull(AffyID)
set02 <- set01[featureNames(set01) %in% genes, ]
We need
to define the model structure for the contrasts we intend to make. Here
we will define 3 major contrast groups; 1) effect of time in placebo
group, 2) effect of time in felzartamab group, and 3) the interactive
effect of time and treatment. We do this for each time window (baseline
to week 24, week 24 to week 52, and baseline to week 52). We also
include a molecular estimate for % cortex as a covariate to adjust for
difference in biopsy composition.
# DEFINE FACTOR FOR CONTRASTS ####
Felzartamab_Followup <- set02$Felzartamab_Followup %>% droplevels()
cortex <- set02$Cortexprob
# DESIGN ####
design <- model.matrix(~ 0 + Felzartamab_Followup + cortex)
# CONTRAST DESIGN week24 - baseline ####
contrast_interaction_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
contrast_felzartamab_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
levels = design
)
contrast_placebo_01 <- makeContrasts(
"x = (Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
# CONTRAST DESIGN week52 - week24 ####
contrast_interaction_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
levels = design
)
contrast_felzartamab_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupWeek24_Felzartamab)/2",
levels = design
)
contrast_placebo_02 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupWeek24_Placebo)/2",
levels = design
)
# CONTRAST DESIGN week52 - baseline####
contrast_interaction_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
contrast_felzartamab_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2",
levels = design
)
contrast_placebo_03 <- makeContrasts(
"x = (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2",
levels = design
)
We can
check the specified contrast matrix for interaction terms. This is for
the baseline to week 24 window.
contrast_interaction_01 %>% as_tibble(rownames = "level")
## # A tibble: 7 x 2
## level x = (Felzartamab_FollowupWeek~1
## <chr> <dbl>
## 1 Felzartamab_FollowupBaseline_Placebo 0.5
## 2 Felzartamab_FollowupWeek24_Placebo -0.5
## 3 Felzartamab_FollowupWeek52_Placebo 0
## 4 Felzartamab_FollowupBaseline_Felzartamab -0.5
## 5 Felzartamab_FollowupWeek24_Felzartamab 0.5
## 6 Felzartamab_FollowupWeek52_Felzartamab 0
## 7 cortex 0
## # i abbreviated name:
## # 1: `x = (Felzartamab_FollowupWeek24_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 -(Felzartamab_FollowupWeek24_Placebo-Felzartamab_FollowupBaseline_Placebo)/2`
and this
is for the baseline to week 52 window.
contrast_interaction_03 %>% as_tibble(rownames = "level")
## # A tibble: 7 x 2
## level x = (Felzartamab_FollowupWeek~1
## <chr> <dbl>
## 1 Felzartamab_FollowupBaseline_Placebo 0.5
## 2 Felzartamab_FollowupWeek24_Placebo 0
## 3 Felzartamab_FollowupWeek52_Placebo -0.5
## 4 Felzartamab_FollowupBaseline_Felzartamab -0.5
## 5 Felzartamab_FollowupWeek24_Felzartamab 0
## 6 Felzartamab_FollowupWeek52_Felzartamab 0.5
## 7 cortex 0
## # i abbreviated name:
## # 1: `x = (Felzartamab_FollowupWeek52_Felzartamab-Felzartamab_FollowupBaseline_Felzartamab)/2 - (Felzartamab_FollowupWeek52_Placebo-Felzartamab_FollowupBaseline_Placebo)/2`
We now
fit the models we’ve specified
# FIT BLOCK week24 - baseline LIMMA MODEL ####
fit_interaction_1 <- limma::lmFit(set02, design)
cfit_interaction_1 <- limma::contrasts.fit(fit_interaction_1, contrast_interaction_01)
ebayes_interaction_1 <- limma::eBayes(cfit_interaction_1)
tab_interaction_1 <- limma::topTable(ebayes_interaction_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_1$stdev.unscaled * sqrt(ebayes_interaction_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_1 <- limma::lmFit(set02, design)
cfit_felzartamab_1 <- limma::contrasts.fit(fit_felzartamab_1, contrast_felzartamab_01)
ebayes_felzartamab_1 <- limma::eBayes(cfit_felzartamab_1)
tab_felzartamab_1 <- limma::topTable(ebayes_felzartamab_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_1$stdev.unscaled * sqrt(ebayes_felzartamab_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_1 <- limma::lmFit(set02, design)
cfit_placebo_1 <- limma::contrasts.fit(fit_placebo_1, contrast_placebo_01)
ebayes_placebo_1 <- limma::eBayes(cfit_placebo_1)
tab_placebo_1 <- limma::topTable(ebayes_placebo_1, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_1$stdev.unscaled * sqrt(ebayes_placebo_1$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
# FIT BLOCK week52 - week24 LIMMA MODEL ####
fit_interaction_2 <- limma::lmFit(set02, design)
cfit_interaction_2 <- limma::contrasts.fit(fit_interaction_2, contrast_interaction_02)
ebayes_interaction_2 <- limma::eBayes(cfit_interaction_2)
tab_interaction_2 <- limma::topTable(ebayes_interaction_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_2$stdev.unscaled * sqrt(ebayes_interaction_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_2 <- limma::lmFit(set02, design)
cfit_felzartamab_2 <- limma::contrasts.fit(fit_felzartamab_2, contrast_felzartamab_02)
ebayes_felzartamab_2 <- limma::eBayes(cfit_felzartamab_2)
tab_felzartamab_2 <- limma::topTable(ebayes_felzartamab_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_2$stdev.unscaled * sqrt(ebayes_felzartamab_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_2 <- limma::lmFit(set02, design)
cfit_placebo_2 <- limma::contrasts.fit(fit_placebo_2, contrast_placebo_02)
ebayes_placebo_2 <- limma::eBayes(cfit_placebo_2)
tab_placebo_2 <- limma::topTable(ebayes_placebo_2, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_2$stdev.unscaled * sqrt(ebayes_placebo_2$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
# FIT BLOCK week52 - baseline LIMMA MODEL ####
fit_interaction_3 <- limma::lmFit(set02, design)
cfit_interaction_3 <- limma::contrasts.fit(fit_interaction_3, contrast_interaction_03)
ebayes_interaction_3 <- limma::eBayes(cfit_interaction_3)
tab_interaction_3 <- limma::topTable(ebayes_interaction_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_interaction_3$stdev.unscaled * sqrt(ebayes_interaction_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_felzartamab_3 <- limma::lmFit(set02, design)
cfit_felzartamab_3 <- limma::contrasts.fit(fit_felzartamab_3, contrast_felzartamab_03)
ebayes_felzartamab_3 <- limma::eBayes(cfit_felzartamab_3)
tab_felzartamab_3 <- limma::topTable(ebayes_felzartamab_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_felzartamab_3$stdev.unscaled * sqrt(ebayes_felzartamab_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
fit_placebo_3 <- limma::lmFit(set02, design)
cfit_placebo_3 <- limma::contrasts.fit(fit_placebo_3, contrast_placebo_03)
ebayes_placebo_3 <- limma::eBayes(cfit_placebo_3)
tab_placebo_3 <- limma::topTable(ebayes_placebo_3, adjust = "fdr", sort.by = "p", number = "all") %>%
as_tibble(rownames = "AffyID") %>%
mutate(se = (ebayes_placebo_3$stdev.unscaled * sqrt(ebayes_placebo_3$s2.post)) %>% as.vector(), .after = logFC) %>%
right_join(affymap219 %>% dplyr::select(AffyID, Symb, Gene, PBT), ., by = "AffyID")
The
results are now organized for tabulation
# MERGE BLOCKS ####
tab_block_1 <- tab_interaction_1 %>%
left_join(
tab_felzartamab_1 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_1 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
tab_block_2 <- tab_interaction_2 %>%
left_join(
tab_felzartamab_2 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_2 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
tab_block_3 <- tab_interaction_3 %>%
left_join(
tab_felzartamab_3 %>%
dplyr::select(AffyID, logFC) %>%
rename(flogFC = logFC),
by = "AffyID"
) %>%
left_join(
tab_placebo_3 %>%
dplyr::select(AffyID, logFC) %>%
rename(plogFC = logFC),
by = "AffyID"
) %>%
relocate(plogFC, flogFC, .before = logFC) %>%
arrange(P.Value)
# CALCULATE MEAN GENE EXPRESSION FOR EACH PROBE BETWEEN GROUPINGS ####
means <- fit_interaction_1 %>%
avearrays() %>%
as_tibble(rownames = "AffyID") %>%
mutate_if(is.numeric, ~ 2^. %>% round(0)) %>%
rename_at(vars(contains("Felz")), ~ str_remove(., "Felzartamab_Followup")) %>%
dplyr::select(-contains("Week12"), -any_of(c("cortex")))
# FORMAT TOPTABLES ####
table_interaction_1 <- tab_block_1 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
table_interaction_2 <- tab_block_2 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
table_interaction_3 <- tab_block_3 %>%
arrange(P.Value) %>%
mutate_at(c("P.Value", "adj.P.Val"), as.numeric) %>%
tibble() %>%
left_join(means, by = "AffyID") %>%
dplyr::select(
AffyID, Symb, Gene, PBT,
t, plogFC, flogFC, logFC, se, P.Value, adj.P.Val,
all_of(colnames(means))
) %>%
dplyr::rename(
p = P.Value,
FDR = adj.P.Val
)
deg <- tibble(
design = c(
"Baseline_vs_Week24",
"Week24_vs_Week52",
"Baseline_vs_Week52"
),
table = list(
table_interaction_1,
table_interaction_2,
table_interaction_3
)
)
# FORMAT TABLES TO MAKE FLEXTABLES ####
table_deg <- deg %>%
expand_grid(direction = c("all", "increased", "decreased")) %>%
relocate(direction, .after = design) %>%
mutate(
gene_tables = pmap(
list(direction, table),
function(direction, table) {
if (direction == "increased") {
table <- table %>%
dplyr::filter(logFC > 0) %>%
arrange(p)
} else if (direction == "decreased") {
table <- table %>%
dplyr::filter(logFC < 0) %>%
arrange(p)
}
table %>%
dplyr::select(-t, -se, -FDR, -contains("AffyID"), -contains("MMDx")) %>%
dplyr::slice(1:20) %>%
mutate(
Gene = Gene %>% str_remove("///.*"),
plogFC = plogFC %>% round(2),
flogFC = flogFC %>% round(2),
logFC = logFC %>% round(2),
p = case_when(
p < 0.0001 ~ p %>% formatC(digits = 0, format = "e"),
TRUE ~ p %>% formatC(digits = 4, format = "f")
)
)
}
)
)
# GLOBAL PARAMETERS FOR FLEXTABLES ####
header1 <- c(
"Gene\nsymbol", "Gene", "PBT",
rep("Differential expression", 4),
rep("Mean expression by group", 6)
)
header2 <- c(
"Gene\nsymbol", "Gene", "PBT",
"\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
"\u394\u394\nlogFC", "\u394\u394\nP",
rep("Placebo", 3), rep("Felzartamab", 3)
)
header3 <- c(
"Gene\nsymbol", "Gene", "PBT",
"\u394\nplacebo\nlogFC", "\u394\nfelzartamab\nlogFC",
"\u394\u394\nlogFC", "\u394\u394\nP",
"Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)", "Baseline\n(N=10)", "Week24\n(N=10)", "Week52\n(N=10)"
)
# MAKE FORMATTED FLEXTABLES ####
flextables <- table_deg %>%
mutate(
flextable = pmap(
list(design, gene_tables),
function(design, gene_tables) {
if (design == "Baseline_vs_Week24") {
title <- paste("Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
} else if (design == "Week24_vs_Week52") {
title <- paste("Table i. Top 20 differentially expressed genes between week24 and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
} else if (design == "Baseline_vs_Week52") {
title <- paste("Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value)", sep = "")
}
gene_tables %>%
arrange(logFC) %>%
flextable::flextable() %>%
flextable::delete_part("header") %>%
flextable::add_header_row(top = TRUE, values = header3) %>%
flextable::add_header_row(top = TRUE, values = header2) %>%
flextable::add_header_row(top = TRUE, values = header1) %>%
flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
flextable::merge_v(part = "header") %>%
flextable::merge_h(part = "header") %>%
flextable::border_remove() %>%
flextable::border(part = "header", border = officer::fp_border()) %>%
flextable::border(part = "body", border = officer::fp_border()) %>%
flextable::align(align = "center") %>%
flextable::align(align = "center", part = "header") %>%
flextable::font(fontname = "Arial", part = "all") %>%
flextable::fontsize(size = 7, part = "all") %>%
flextable::fontsize(size = 7, part = "footer") %>%
flextable::fontsize(i = 1, size = 12, part = "header") %>%
flextable::bold(part = "header") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::padding(padding = 0, part = "all") %>%
flextable::autofit()
}
)
)
Genes
suppressed by felzartamab in week 24 compared to baseline
(Note
that all ΔΔFDR are > 0.05. This is a consequence of differential
expression across a large selection of genes (5,267 in this case) with
small samples sizes (i.e., n = 10 by treatment group). Thus, treatment
effects on individual genes should not be interpreted. Instead, geneset
enrichment analysis (GSEA) will be used to assess bulk
trends).
flextables %>%
dplyr::filter(design == "Baseline_vs_Week24", direction == "decreased") %>%
pull(flextable) %>%
pluck(1)
Table i. Top 20 differentially expressed genes between baseline and week24 in biopsies from placebo and Felzartamab treated patients (by P-value) | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Gene | Gene | PBT | Differential expression | Mean expression by group | ||||||||
Δ | Δ | ΔΔ | ΔΔ | Placebo | Felzartamab | |||||||
Baseline | Week24 | Week52 | Baseline | Week24 | Week52 | |||||||
CXCL11 | chemokine (C-X-C motif) ligand 11 | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.16 | -1.01 | -0.85 | 0.0176 | 494 | 395 | 372 | 617 | 153 | 317 |
CXCL9 | chemokine (C-X-C motif) ligand 9 | ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT | -0.14 | -0.94 | -0.79 | 0.0290 | 1,344 | 1,102 | 1,384 | 1,684 | 459 | 1,036 |
CXCL10 | chemokine (C-X-C motif) ligand 10 | ABMR-RAT,GRIT1,GRIT3,RAT,Rej-RAT | -0.21 | -0.98 | -0.77 | 0.0228 | 844 | 632 | 720 | 1,063 | 274 | 615 |
IDO1 | indoleamine 2,3-dioxygenase 1 | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.09 | -0.77 | -0.68 | 0.0103 | 338 | 298 | 296 | 409 | 140 | 251 |
LYZ | lysozyme | IRITD5,QCMAT | -0.01 | -0.67 | -0.66 | 0.0386 | 1,676 | 1,658 | 1,884 | 2,244 | 884 | 1,184 |
FCGR3A | Fc fragment of IgG, low affinity IIIa, receptor (CD16a) | ABMR-RAT,GRIT2,IRRAT950,RAT,Rej-RAT | 0.06 | -0.58 | -0.64 | 0.0167 | 550 | 602 | 613 | 571 | 256 | 352 |
GZMB | granzyme B | ABMR-RAT,QCAT,RAT,Rej-RAT | -0.01 | -0.62 | -0.61 | 0.0019 | 149 | 147 | 130 | 143 | 61 | 114 |
GBP1 | guanylate binding protein 1, interferon-inducible | ABMR-RAT,GRIT3,RAT,Rej-RAT | -0.03 | -0.61 | -0.59 | 0.0096 | 678 | 655 | 643 | 909 | 389 | 556 |
GNLY | granulysin | ABMR-RAT,DSAST,QCAT,RAT,Rej-RAT | -0.18 | -0.75 | -0.57 | 0.0285 | 192 | 149 | 134 | 165 | 58 | 100 |
SH2D1B | SH2 domain containing 1B | ABMR-RAT,DSAST,NKB,RAT | 0.07 | -0.44 | -0.51 | 0.0010 | 38 | 42 | 34 | 40 | 22 | 33 |
FCGR1A | Fc fragment of IgG, high affinity Ia, receptor (CD64) | GRIT2,IRRAT950,RAT,TCMR-RAT | 0.09 | -0.41 | -0.50 | 0.0304 | 131 | 148 | 161 | 152 | 86 | 86 |
KLRD1 | killer cell lectin-like receptor subfamily D, member 1 | ABMR-RAT,RAT,Rej-RAT | -0.06 | -0.49 | -0.43 | 0.0295 | 110 | 101 | 96 | 126 | 63 | 102 |
EPB41L3 | erythrocyte membrane protein band 4.1-like 3 | 0.22 | -0.20 | -0.42 | 0.0388 | 209 | 286 | 258 | 339 | 258 | 283 | |
LOC339059 | uncharacterized LOC339059 | 0.26 | -0.15 | -0.41 | 0.0203 | 51 | 73 | 55 | 66 | 53 | 59 | |
APOL1 | apolipoprotein L1 | ABMR-RAT,GRIT3,RAT,Rej-RAT | 0.07 | -0.32 | -0.39 | 0.0205 | 626 | 687 | 744 | 754 | 483 | 636 |
LYPD5 | LY6/PLAUR domain containing 5 | ABMR-RAT,RAT,Rej-RAT | 0.04 | -0.33 | -0.37 | 0.0255 | 27 | 28 | 27 | 30 | 19 | 26 |
TRAV12-2 | T cell receptor alpha variable 12-2 | 0.25 | -0.10 | -0.35 | 0.0173 | 39 | 55 | 39 | 49 | 43 | 43 | |
IL18RAP | interleukin 18 receptor accessory protein | 0.02 | -0.32 | -0.34 | 0.0233 | 38 | 39 | 41 | 42 | 27 | 33 | |
PRF1 | perforin 1 (pore forming protein) | ABMR-RAT,QCAT,RAT,Rej-RAT | -0.02 | -0.35 | -0.33 | 0.0382 | 152 | 148 | 137 | 145 | 89 | 123 |
PMM2 | phosphomannomutase 2 | 0.20 | -0.07 | -0.27 | 0.0380 | 38 | 50 | 41 | 46 | 42 | 49 | |
Genes
suppressed by felzartamab in week 52 compared to baseline
(Note
that all ΔΔFDR are > 0.05. This is a consequence of differential
expression across a large selection of genes (5,267 in this case) with
small samples sizes (i.e., n = 10 by treatment group). Thus, treatment
effects on individual genes should not be interpreted. Instead, geneset
enrichment analysis (GSEA) will be used to assess bulk
trends).
flextables %>%
dplyr::filter(design == "Baseline_vs_Week52", direction == "decreased") %>%
pull(flextable) %>%
pluck(1)
Table i. Top 20 differentially expressed genes between baseline and week52 in biopsies from placebo and Felzartamab treated patients (by P-value) | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Gene | Gene | PBT | Differential expression | Mean expression by group | ||||||||
Δ | Δ | ΔΔ | ΔΔ | Placebo | Felzartamab | |||||||
Baseline | Week24 | Week52 | Baseline | Week24 | Week52 | |||||||
HNRNPK | heterogeneous nuclear ribonucleoprotein K | 0.28 | -0.38 | -0.66 | 0.0055 | 818 | 1,044 | 1,203 | 1,096 | 1,059 | 649 | |
COL1A1 | collagen, type I, alpha 1 | COLA1356,FICOL,IRITD5,LivGST_UP | 0.38 | -0.26 | -0.64 | 0.0075 | 310 | 333 | 526 | 329 | 320 | 230 |
NFIX | nuclear factor I/X (CCAAT-binding transcription factor) | CT1 | 0.42 | -0.15 | -0.57 | 0.0012 | 162 | 179 | 289 | 221 | 167 | 180 |
RHOA | ras homolog family member A | IRRAT950 | 0.28 | -0.29 | -0.57 | 0.0049 | 691 | 803 | 1,017 | 923 | 809 | 620 |
COL1A2 | collagen, type I, alpha 2 | FICOL,IRITD5,LivGST_UP | 0.26 | -0.31 | -0.57 | 0.0118 | 1,238 | 1,152 | 1,769 | 1,534 | 1,194 | 997 |
ITGB3 | integrin beta 3 | IRRAT30,IRRAT950 | 0.34 | -0.21 | -0.55 | 0.0100 | 105 | 137 | 169 | 110 | 103 | 82 |
LOC100129518 | uncharacterized LOC100129518 | GRIT3,IRRAT950 | 0.25 | -0.22 | -0.47 | 0.0050 | 429 | 546 | 605 | 527 | 461 | 390 |
ACTB | actin, beta | IRRAT950 | 0.15 | -0.31 | -0.47 | 0.0077 | 6,451 | 7,352 | 7,993 | 7,274 | 5,662 | 4,720 |
SPARC | secreted protein, acidic, cysteine-rich (osteonectin) | IRITD3 | 0.16 | -0.30 | -0.46 | 0.0100 | 2,411 | 2,532 | 3,011 | 2,343 | 2,094 | 1,550 |
DNAJA4 | DnaJ (Hsp40) homolog, subfamily A, member 4 | IRRAT950 | 0.27 | -0.14 | -0.42 | 0.0095 | 49 | 68 | 72 | 68 | 61 | 56 |
SLFN5 | schlafen family member 5 | 0.19 | -0.20 | -0.40 | 0.0050 | 39 | 44 | 51 | 49 | 38 | 37 | |
CYR61 | cysteine-rich, angiogenic inducer, 61 | IRITD3,IRRAT950 | 0.14 | -0.26 | -0.40 | 0.0074 | 271 | 268 | 330 | 338 | 234 | 236 |
SDCBP | syndecan binding protein | 0.24 | -0.14 | -0.38 | 0.0120 | 1,135 | 1,404 | 1,582 | 1,478 | 1,358 | 1,224 | |
SBSPON | somatomedin B and thrombospondin type 1 domain containing | 0.31 | -0.08 | -0.38 | 0.0123 | 217 | 250 | 332 | 188 | 249 | 169 | |
PAFAH1B2 | platelet-activating factor acetylhydrolase 1b, catalytic subunit 2 (30kDa) | 0.23 | -0.13 | -0.36 | 0.0099 | 309 | 372 | 423 | 386 | 396 | 323 | |
TMEM2 | transmembrane protein 2 | 0.10 | -0.26 | -0.36 | 0.0118 | 270 | 307 | 311 | 318 | 259 | 222 | |
CALM2 | calmodulin 2 (phosphorylase kinase, delta) | CT1,IRITD3 | 0.12 | -0.25 | -0.36 | 0.0138 | 1,304 | 1,427 | 1,538 | 1,623 | 1,342 | 1,155 |
WNK1 | WNK lysine deficient protein kinase 1 | 0.26 | -0.08 | -0.34 | 0.0104 | 703 | 884 | 1,012 | 894 | 996 | 799 | |
CEBPB | CCAAT/enhancer binding protein (C/EBP), beta | IRITD3 | 0.11 | -0.22 | -0.33 | 0.0123 | 419 | 452 | 488 | 417 | 333 | 307 |
STAT6 | signal transducer and activator of transcription 6, interleukin-4 induced | 0.09 | -0.20 | -0.30 | 0.0144 | 671 | 710 | 765 | 746 | 616 | 563 | |
Geneset
enrichment analysis (GSEA) is applied to identify pathways affected by
felzartamab treatment. Multiple libraries were screened to encompass a
broad spectrum of biological pathways, including Gene Ontology (GO),
Disease Ontology Semantic and Enrichment analysis (DOSE), Kyoto
Encyclopedia of Genes and Genomes (KEGG), Reactome, Molecular Signatures
Database (MsigDB), and WikiPathways (wiki) libraries. An additional
de-novo injury geneset was also assesed. This injury geneset was based
on single-nuclei transcriptomics assessed from healthy and injury kidney
allografts. Geneset enrichment was chosen specifcally because it is a
rank-based test and does not require cutoffs based on p-values or
t-values. Nonetheless, here we decided apriori that we would pre-screen
genes based on their p-values from differential expression analyses
described in Section 3. Note that there are many ways for carrying out
pathway and functional enrichment analyses. Also note that gene
libraries are frequently updated and you may find slightly different
results based on the packge version being used. For example, here we
encountered such a difference in the DOSE library, where versions later
than 3.30.1 identified two different pathways. Thus, we have included a
tarball of the package version used for the original analyses in this
paper for convenient installation and replication of the presented
findings.
References:
• Yu, G., Wang, L.G., Han, Y. & He, Q.Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284-287 (2012).
• Wu, T., et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2, 100141 (2021).
• Hinze, C., et al. Single-cell transcriptomics reveals common epithelial response patterns in human acute kidney injury. Genome Med 14, 103 (2022).
First we
need to load the necessary R packages and data.
# HOUSEKEEPING ####
# CRAN libraries
library(tidyverse)
library(flextable)
library(officer)
library(msigdbr)
# Bioconductor libraries
# DOSE v 3.30.1 required for reproducibility
# install.packages("natmed/data/DOSE_3.30.1.tar.gz", repos = NULL, type = "source")
library(DOSE)
library(ReactomePA)
library(Biobase)
library(clusterProfiler)
library(org.Hs.eg.db)
library(pRoloc)
## Warning in fun(libname, pkgname): mzR has been built against a different Rcpp version (1.0.12)
## than is installed on your system (1.0.13). This might lead to errors
## when loading mzR. If you encounter such issues, please send a report,
## including the output of sessionInfo() to the Bioc support forum at
## https://support.bioconductor.org/. For details see also
## https://github.com/sneumann/mzR/wiki/mzR-Rcpp-compiler-linker-issue.
# Custom operators and functions
"%nin%" <- function(a, b) match(a, b, nomatch = 0) == 0
# load DE data
load("C:/R/CD38-effect-of-treatment/natmed/results/deg.RData")
# Seed for reproducibility
set.seed(42)
Geneset
enrichment analyses using the ClusterProfiler package have simlar, but
not identical implementation. This section will carry out GSEA on all
described libraries, after which the results will be collated and
summarized.
GO
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_go <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::gseGO(
gene = genes_gsea, ont = "BP", OrgDb = org.Hs.eg.db,
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# SIMPLIFY GO TERMS ####
set.seed(42)
gsea_go_simplified <- gsea_go %>%
mutate(gsea_simplified = map(gsea, clusterProfiler::simplify))
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_go_simplified_formatted <- gsea_go_simplified %>%
mutate(
gsea_tables = map(
gsea_simplified,
function(gsea_simplified) {
gsea_simplified %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_go_tables <- gsea_go_simplified_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_go <- gsea_go_tables %>%
mutate(db = "GO", .before = 1)
names(gsea_go$gsea_flextables) <- gsea_go$design
DOSE
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = clusterProfiler::bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_dose <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
DOSE::gseDO(
gene = genes_gsea,
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_dose_formatted <- gsea_dose %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"))
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_dose_tables <- gsea_dose_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_dose <- gsea_dose_tables %>%
mutate(db = "DOSE", .before = 1)
names(gsea_dose$gsea_flextables) <- gsea_dose$design
KEGG
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_kegg <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::gseKEGG(
gene = genes_gsea, organism = "hsa",
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_kegg_formatted <- gsea_kegg %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_kegg_tables <- gsea_kegg_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_kegg <- gsea_kegg_tables %>%
mutate(db = "kegg", .before = 1)
names(gsea_kegg$gsea_flextables) <- gsea_kegg$design
MSigDB
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# DEFINE MSigDB PATHWAYS ####
gene_sets <- msigdbr(category = "H")
map <- gene_sets[, c("gs_name", "entrez_gene")]
map$entrez_gene <- as.character(map$entrez_gene)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_msigdb <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::GSEA(
gene = genes_gsea, TERM2GENE = map,
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_msigdb_formatted <- gsea_msigdb %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_msigdb_tables <- gsea_msigdb_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
# slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_msigdb <- gsea_msigdb_tables %>%
mutate(db = "msigdb", .before = 1)
names(gsea_msigdb$gsea_flextables) <- gsea_msigdb$design
REACTOME
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_reactome <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
ReactomePA::gsePathway(
gene = genes_gsea, organism = "human",
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_reactome_formatted <- gsea_reactome %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_reactome_tables <- gsea_reactome_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_reactome <- gsea_reactome_tables %>%
mutate(db = "Reactome", .before = 1)
names(gsea_reactome$gsea_flextables) <- gsea_reactome$design
wikiPathways
library
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_wiki <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::gseWP(
gene = genes_gsea, organism = "Homo sapiens",
minGSSize = 10, maxGSSize = 200,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_wiki_formatted <- gsea_wiki %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = factor(Description, levels = Description, ordered = TRUE)
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_wiki_tables <- gsea_wiki_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
# slice_min(pvalue, n = 20) %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_wiki <- gsea_wiki_tables %>%
mutate(db = "wiki", .before = 1)
names(gsea_wiki$gsea_flextables) <- gsea_wiki$design
AKI
geneset
# load gene lists
load("C:/R/CD38-effect-of-treatment/natmed/data/hinze2022_injury_markers.RData")
# DEFINE INJURY PATHWAYS ####
map_aki <- injury_markers %>% dplyr::select(gs_name, entrez_gene)
map_aki_joined <- injury_markers %>%
dplyr::select(gs_name, entrez_gene) %>%
mutate(gs_name = "AKI")
# WRANGLE GENE DATA FOR GSEA ####
data <- deg %>%
mutate(
genes_gsea = map(
table,
function(table) {
table %>%
dplyr::filter(p < 0.05) %>%
arrange(t %>% dplyr::desc()) %>%
dplyr::mutate(
ENTREZID = bitr(Symb,
fromType = "SYMBOL",
toType = c("ENTREZID"),
OrgDb = org.Hs.eg.db,
drop = FALSE
) %>% pull(ENTREZID)
) %>%
drop_na(ENTREZID) %>%
dplyr::select(t, ENTREZID) %>%
pull(t, ENTREZID)
}
)
)
# GENESET ENRICHMENT ANALYSES (GSEA) ####
set.seed(42)
gsea_aki <- data %>%
mutate(
gsea = map(
genes_gsea,
function(genes_gsea) {
clusterProfiler::GSEA(
gene = genes_gsea, TERM2GENE = map_aki_joined,
minGSSize = 5, maxGSSize = Inf,
pvalueCutoff = 0.001, pAdjustMethod = "fdr", seed = TRUE
) %>% clusterProfiler::setReadable(OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_aki_formatted <- gsea_aki %>%
mutate(
gsea_tables = map(
gsea,
function(gsea) {
gsea %>%
as_tibble() %>%
arrange(pvalue) %>%
mutate(
sign = case_when(NES < 0 ~ "Down Regulated", NES > 0 ~ "Up Regulated"),
Description = "Cellular response to acute kidney injury"
)
}
)
)
# FORMAT TABLES FOR GENESET ENRICHMENT ANALYSES (GSEA) ####
gsea_aki_tables <- gsea_aki_formatted %>%
mutate(
gsea_flextables = map(
gsea_tables,
function(gsea_tables) {
gsea_tables %>%
mutate(
NES = NES %>% round(2),
pvalue = case_when(
pvalue < 0.0001 ~ pvalue %>% formatC(digits = 0, format = "e"),
TRUE ~ pvalue %>% formatC(digits = 4, format = "f")
),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(Description, setSize, NES, pvalue, FDR, dplyr::any_of("core_enrichment")) %>%
flextable::flextable() %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core_enrichment", align = "left", part = "body") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::autofit()
}
)
)
# PREPARE THE RESULTS FOR EXPORT ####
gsea_aki <- gsea_aki_tables %>%
mutate(db = "Hinze", .before = 1)
names(gsea_aki$gsea_flextables) <- gsea_aki$design
GSEA
results from each library can now be combined and summarised. Note that
pathway annotations need to be interpreted in light of their
developement in combination with how GSEA determines their enrichment,
and the relevance of enriched pathways must be evaluated in the context
of transplantation. For our interpretation, we avail to the pathway
annotation as well as core enrichment genes contributing to enrichment
via GSEA. We have simplified the interpretation of enriched pathways
based on the overlap in biological function of their core enrichement
genes.
# JOIN THE GSEA RESULTS ####
gsea <- purrr::reduce(
list(
gsea_go,
gsea_msigdb,
gsea_wiki,
gsea_dose,
gsea_kegg,
gsea_reactome,
gsea_aki
), bind_rows
) %>%
mutate(keep = map_dbl(gsea_tables, nrow)) %>%
dplyr::filter(keep > 0) %>%
dplyr::select(db, design, genes_gsea, gsea_tables) %>%
mutate(gsea_tables = map(gsea_tables, mutate, Description = Description %>% as.character()))
# DEFINE INTERPREATION BASED ON PATHWAYS ####
pathway_interpretation <- gsea %>%
dplyr::select(design, gsea_tables) %>%
unnest(everything()) %>%
mutate(
interpretation = case_when(
ID %in% c(
"GO:0006952",
"GO:0009605",
"GO:0009607",
"GO:0043207",
"GO:0044419",
"GO:0051707",
"GO:0006955",
"GO:0002376",
"DOID:2914",
"DOID:77",
"DOID:1579",
"R-HSA-168256"
) & design == "Baseline_vs_Week24" ~ "immune response",
ID %in% c(
"R-HSA-168256",
"R-HSA-168249",
"R-HSA-1643685",
"R-HSA-5663205"
) & design == "Baseline_vs_Week52" ~ "immune-related response to injury",
ID %in% c("AKI", "R-HSA-109582") ~ "response to injury",
ID %in% c(
"R-HSA-9006934",
"R-HSA-597592",
"R-HSA-5653656",
"R-HSA-392499"
) ~ "response to injury",
ID %in% c(
"R-HSA-2262752",
"R-HSA-8953897"
) ~ "response to injury",
)
) %>%
dplyr::select(design, ID, Description, interpretation)
# WRANGLE THE GSEA TABLES ####
gsea_tables <- gsea %>%
mutate(
gsea_tables = pmap(
list(design, gsea_tables),
function(design, gsea_tables) {
gsea_tables %>%
mutate(
design = design,
Description = Description %>% as.character()
) %>%
left_join(pathway_interpretation, by = c("ID", "Description", "design")) %>%
dplyr::select(ID, Description, interpretation, setSize, NES, p.adjust, core_enrichment) %>%
distinct(ID, setSize, interpretation, .keep_all = TRUE)
}
)
) %>%
unnest(gsea_tables) %>%
arrange(p.adjust) %>%
dplyr::filter(p.adjust < 0.001) %>%
nest(.by = design, gsea_table = c(-design, -genes_gsea)) %>%
left_join(gsea %>% dplyr::select(design, genes_gsea) %>% distinct(design, .keep_all = TRUE),
by = "design"
)
# MAKE FLEXTABLES ####
gsea_flextables <- gsea_tables %>%
mutate(
flextable = pmap(
list(design, gsea_table),
function(design, gsea_table) {
cellWidths <- c(2, 3, 7, 4, 1.5, 1.5, 1.5, 13) # for individual tables up or down
if (design == "Baseline_vs_Week24") {
title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week24 (by FDR)", sep = "")
} else if (design == "Week24_vs_Week52") {
title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between week24 and week52 (by FDR)", sep = "")
} else if (design == "Baseline_vs_Week52") {
title <- paste("Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week52 (by FDR)", sep = "")
}
gsea_table %>%
mutate(
library = db,
pathway = Description,
"core enrichment genes" = core_enrichment %>% str_replace_all("/", ", "),
"n genes" = setSize,
NES = NES %>% round(2),
FDR = case_when(
p.adjust < 0.0001 ~ p.adjust %>% formatC(digits = 0, format = "e"),
TRUE ~ p.adjust %>% formatC(digits = 4, format = "f")
),
) %>%
dplyr::select(library, ID, pathway, interpretation, "n genes", NES, FDR, "core enrichment genes") %>%
flextable::flextable() %>%
flextable::add_header_row(top = TRUE, values = rep(title, ncol_keys(.))) %>%
flextable::merge_v(part = "header") %>%
flextable::merge_h(part = "header") %>%
flextable::border_remove() %>%
flextable::border(border = officer::fp_border(), part = "all") %>%
flextable::padding(padding = 0) %>%
flextable::padding(j = "core enrichment genes", padding.left = 5) %>%
flextable::align(align = "center", part = "all") %>%
flextable::align(j = "core enrichment genes", align = "left", part = "body") %>%
flextable::font(fontname = "Arial", part = "all") %>%
flextable::fontsize(size = 10, part = "body") %>%
flextable::fontsize(size = 12, part = "header") %>%
flextable::bg(bg = "white", part = "all") %>%
flextable::width(width = cellWidths, unit = "cm")
}
)
)
Summary
of all enriched pathways at week 24 (compared to
baseline)
# PRINT FLETABLES TO POWERPOINT ####
gsea_flextables %>%
dplyr::filter(design == "Baseline_vs_Week24") %>%
pull(flextable) %>%
pluck(1)
Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week24 (by FDR) | |||||||
|---|---|---|---|---|---|---|---|
library | ID | pathway | interpretation | n genes | NES | FDR | core enrichment genes |
GO | GO:0006952 | defense response | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0009605 | response to external stimulus | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0009607 | response to biotic stimulus | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0043207 | response to external biotic stimulus | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0044419 | biological process involved in interspecies interaction between organisms | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0051707 | response to other organism | immune response | 20 | -3.07 | 3e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0006955 | immune response | immune response | 23 | -2.93 | 5e-05 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, TRAV12-2, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
GO | GO:0002376 | immune system process | immune response | 24 | -2.75 | 0.0002 | PTPRC, TRIM4, LYZ, PRF1, FCGR1A, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, APOL1, CXCL11, TRAV12-2, FCGR3A, IDO1, GBP1, GZMB, SH2D1B |
Reactome | R-HSA-168256 | Immune System | immune response | 14 | -2.51 | 0.0005 | PIK3AP1, PTPRC, TRIM4, KLRF1, LYZ, FCGR1A, KLRD1, GNLY, IL18RAP, CXCL10, FCGR3A, GBP1, SH2D1B |
DOSE | DOID:2914 | immune system disease | immune response | 10 | -2.56 | 0.0007 | PRF1, FCGR1A, CXCL9, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB |
DOSE | DOID:77 | gastrointestinal system disease | immune response | 15 | -2.59 | 0.0007 | PTPRC, EPB41L3, LYZ, PRF1, KLRD1, CXCL9, GNLY, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB |
DOSE | DOID:1579 | respiratory system disease | immune response | 10 | -2.56 | 0.0007 | PRF1, FCGR1A, CXCL9, IL18RAP, CXCL10, CXCL11, FCGR3A, IDO1, GZMB |
Summary
of all enriched pathways at week 52 (compared to
baseline)
gsea_flextables %>%
dplyr::filter(design == "Baseline_vs_Week52") %>%
pull(flextable) %>%
pluck(1)
Table Si. Functional enrichment analysis of genes affects by felzartamab treatment between baseline and week52 (by FDR) | |||||||
|---|---|---|---|---|---|---|---|
library | ID | pathway | interpretation | n genes | NES | FDR | core enrichment genes |
Reactome | R-HSA-168249 | Innate Immune System | immune-related response to injury | 20 | -3.18 | 6e-06 | LRRFIP1, C1QB, JUN, ATP6V1A, NCF2, CFI, EP300, FCGR1A, CAPZA1, CAPZA2, ADAM10, HSPA1A, STAT6, CALM2, SDCBP, PAFAH1B2, ACTB, RHOA |
Reactome | R-HSA-168256 | Immune System | immune-related response to injury | 35 | -3.08 | 2e-05 | LRRFIP1, C1QB, CANX, JUN, ATP6V1A, RNF213, NCF2, PIK3AP1, VIM, CFI, EP300, FCGR1A, CAPZA1, CAPZA2, ADAM10, HSPA1A, IFNGR1, STAT6, CALM2, SDCBP, COL1A2, PAFAH1B2, ACTB, COL1A1, RHOA |
Reactome | R-HSA-1643685 | Disease | immune-related response to injury | 32 | -2.92 | 2e-05 | LRRFIP1, PPFIBP1, CANX, JUN, RNF213, CD163, PIK3AP1, TPM3, SIN3A, EP300, FCGR1A, RAN, ADAM10, HSPA1A, IFNGR1, CALM2, CEBPB, SBSPON, COL1A2, ITGB3, ACTB, COL1A1, HNRNPK |
Hinze | AKI | Cellular response to acute kidney injury | response to injury | 39 | -2.75 | 4e-05 | SPTBN1, NRP1, KLF6, HLA-DRB1, DOCK11, DUSP1, CHSY1, LDHA, CDC42SE2, TPM4, PTBP3, ARHGAP29, SYNE2, SERP1, CORO1C, LRRFIP1, PPFIBP1, CANX, JUN, SCOC, RNF213, SPP1, PIK3AP1, VIM, SVIL, TPM3, RAN, ANXA5, ADAM10, HSPA1A, CEBPB, SPARC, ACTB, HNRNPK, SLFN5, RHOA |
Reactome | R-HSA-9006934 | Signaling by Receptor Tyrosine Kinases | response to injury | 19 | -2.59 | 0.0002 | ATP6V1A, SPP1, NCF2, SIN3A, EP300, ADAM10, STAT6, CALM2, COL1A2, SPARC, ITGB3, ACTB, COL1A1, RHOA |
Reactome | R-HSA-597592 | Post-translational protein modification | response to injury | 17 | -2.68 | 0.0002 | CANX, SPP1, PRMT3, PPP6R3, SIN3A, EP300, CAPZA1, CAPZA2, ADAM10, CALM2, SBSPON, ACTB, HNRNPK, RHOA |
Reactome | R-HSA-109582 | Hemostasis | response to injury | 17 | -2.54 | 0.0005 | SIN3A, CAPZA1, CAPZA2, ANXA5, CALM2, COL1A2, SPARC, ITGB3, ACTB, COL1A1, RHOA |
Reactome | R-HSA-392499 | Metabolism of proteins | response to injury | 19 | -2.42 | 0.0005 | PRMT3, PPP6R3, SIN3A, EP300, CAPZA1, CAPZA2, ADAM10, CALM2, SBSPON, ACTB, HNRNPK, RHOA |
Reactome | R-HSA-5653656 | Vesicle-mediated transport | response to injury | 15 | -2.47 | 0.0008 | SCOC, AP4B1, CD163, PPP6R3, CAPZA1, CAPZA2, CALM2, COL1A2, SPARC, PAFAH1B2, ACTB, COL1A1 |